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Abstract 

A  numerical  simulation  has  been  performed  of  the  disturbed  flow  through  and  over  a  two- 
dimensional  array  of  rectangular  buildings  immersed  in  a  neutrally  stratified  deep  rough- 
walled  turbulent  boundary-layer  flow.  The  model  used  for  the  simulation  was  the  steady-state 
Reynolds-averaged  Navier-Stokes  equations  with  linear  and  non-linear  eddy  viscosity 
formulations  for  the  Reynolds  stresses.  The  eddy  viscosity  was  determined  using  a  high- 
Reynolds  number  form  of  the  k-e.  turbulence-closure  model  with  the  boundary  conditions  at 
the  wall  obtained  with  a  standard  wall-function  approach.  The  resulting  system  of  partial 
differential  equations  was  solved  using  the  SIMPLE  algorithm  in  conjunction  with  a  non- 
orthogonal,  colocated,  cell-centered,  finite  volume  procedure.  The  predictive  capabilities  of  the 
high-resolution  computational  fluid  dynamics  (CFD)  simulations  of  urban  flow  are  validated 
against  a  very  detailed  and  comprehensive  wind  tunnel  data  set.  Vertical  profiles  of  the  mean 
streamwise  velocity  and  the  turbulence  kinetic  energy  are  presented  and  compared  to  those 
measured  in  the  wind  tunnel  simulation. 

It  is  found  that  the  performance  of  all  the  turbulence  models  investigated  is  generally 
good — most  of  the  qualitative  features  in  the  disturbed  turbulent  flow  field  through  and  over 
the  building  array  are  correctly  reproduced.  The  quantitative  agreement  is  also  fairly  good 
(especially  for  the  mean  velocity  field).  Overall,  the  non-linear  k-e  model  gave  the  best 
performance  among  four  different  turbulence  closure  models  examined.  The  turbulence  energy 
levels  within  the  street  canyons  and  in  the  exit  region  downstream  of  the  last  building  were 
underestimated  by  all  four  turbulence  closure  models.  This  appears  to  contradict  the 
‘stagnation  point  anomaly’  associated  with  the  standard  k-e  model  which  is  a  result  of  the 
excessive  turbulence  energy  production  due  to  normal  straining.  A  possible  explanation  for 
this  is  the  inability  of  the  present  models  to  account  properly  for  the  effects  of  secondary 
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strains  on  the  turbulence  and/or  for  the  effects  of  large-scale  flapping  of  the  strong  shear  layer 
at  the  canopy  top. 

The  results  of  the  high-resolution  CFD  simulations  have  been  used  to  diagnose  values  of  the 
drag  coefficient  to  be  used  in  a  distributed  drag  force  representation  of  the  obstacles  in  the 
array.  Comparisons  of  the  measured  spatially-averaged  time-mean  mean  velocity  and 
turbulence  kinetic  energy  in  the  array  with  predictions  of  the  disturbed  flow  using  the 
distributed  drag  force  approach  have  been  made. 

©  2003  Elsevier  Ltd.  All  rights  reserved. 
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1.  Introduction 

In  recent  years,  a  great  effort  has  been  made  to  understand  the  impact  of  large 
numbers  of  discrete  bluff  obstacles  (e.g.,  buildings,  trees,  and  other  obstructions  in  a 
prescribed  domain  within  a  large  city)  on  the  atmospheric  boundary-layer  flow 
because  of  its  importance  in  many  aspects  of  meteorology,  wind  engineering,  and 
environmental  science.  In  particular,  how  the  speed,  direction,  and  turbulence  of  the 
wind  changes  as  it  blows  through  a  cluster  of  buildings  must  be  known  if  the  mean 
transport  and  turbulent  diffusion  of  airborne  contaminants  in  its  vicinity  are  to  be 
determined.  Furthermore,  knowledge  of  the  velocity  and  pressure  distributions  near 
the  inlets  and  outlets  of  heating  and  ventilating  systems  of  buildings  is  required  to 
determine  the  egress  of  contaminants  into  buildings.  Indeed,  in  response  to 
increasing  concerns  on  the  use  of  chemical,  biological,  and  radiological  materials 
against  populations  in  urban  areas,  there  has  been  a  concerted  effort  by  a  number  of 
government,  military,  and  civilian  research  laboratories  and  agencies  in  recent  years 
to  develop  models  that  can  predict  the  turbulent  flow  and  dispersion  in  urban  (built- 
up)  areas. 

As  far  as  numerical  modeling  of  flows  round  obstacles  is  concerned,  most  of  the 
effort  has  been  concentrated  on  modeling  the  flow  around  isolated  structures  or 
single  buildings  rather  than  around  groups  of  obstacles  or  building  clusters. 
Prognostic  models  that  solve  the  Reynolds-averaged  Navier  Stokes  (RANS) 
equations  have  been  used  to  describe  the  flow  field  around  an  isolated  cuboid  or 
building  (see  Refs.  [1-4])  or  the  interaction  of  two  buildings  [5].  These  models  used 
the  simple  k—e  turbulence  model  in  the  RANS  equations  to  determine  the  mean  flow 
and  turbulence  and,  hence,  were  basically  non-linear  flow  models  solved  using  either 
finite  difference,  finite  element,  or  finite  volume  solution  techniques.  Computational 
fluid  dynamics  (CFD)  models  have  also  been  applied  to  describe  the  flow  within  a 
specific  urban  structure  (e.g.,  urban  street  canyon).  In  particular,  Hunter  et  al.  [6] 
used  the  k-e  model  of  turbulence  to  describe  the  three-dimensional  characteristics  of 
flow  regimes  within  an  urban  canyon. 

The  most  complicated  flow  problem  is  arguably  that  of  a  large  building  cluster 
consisting  of  a  group  of  buildings  of  roughly  comparable  size  with  the  spacing 
between  buildings  sufficiently  close  that  the  flow  round  any  individual  building  in  the 
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cluster  is  being  influenced  by  all  other  buildings  in  its  neighborhood  (e.g.,  see  the 
discussion  provided  by  Hosker  [7]).  Generally  speaking,  urban  flow  and  dispersion 
within  various  arrays  of  obstacles  are  addressed  on  a  case-by-case  basis.  From  a 
theoretical  point  of  view,  the  problem  appears  to  be  complex  and  unwieldy  and  there 
appears  to  be  no  idealized  starting  point  for  analysis.  From  an  application  point  of 
view,  it  is  easy  to  get  lost  in  the  complexities  of  particular  cases  without  ever  having 
to  see  that  there  might  be  a  common  physical  picture  of  mean  flow  and  turbulence 
characteristics  in  various  building  cluster  configurations.  However,  from  a  numerical 
viewpoint,  it  appears  that  numerical  simulation  models  based  on  RANS  equations 
are  general  enough  to  be  applied  mutatis  mutandi  to  the  problem  of  disturbed  flow 
within  general  building  cluster  configurations  where  mutual  flow  interferences 
between  buildings  need  to  be  taken  into  account.  The  principal  problem  with  the 
application  of  such  CFD  models  to  building  arrays  is  that  there  is  currently  very  little 
benchmark  mean  flow  and  turbulence  data  that  could  be  used  to  compare  with  the 
model  predictions.  It  is  crucial  that  before  any  CFD  model  with  its  attendant 
turbulence  closure  is  applied  to  a  particular  class  of  problems  (e.g.,  flow  through 
building  array),  that  the  model  predictions  are  subjected  first  to  critical  and 
comparative  judgements  against  high-quality  data  sets  before  such  model  predictions 
are  taken  as  “ground  truth”. 

In  this  paper,  we  use  the  k~e  model  of  turbulence  within  the  RANS  modeling 
framework  to  simulate  the  developing  flow  through  and  over  a  two-dimensional 
(2D)  array  of  surface  mounted  obstacles  immersed  within  a  deep  rough-walled 
turbulent  boundary-layer  under  neutrally  stratified  flow.  We  show  a  detailed 
comparison  with  a  wind  tunnel  experiment  of  the  mean  flow  field  and  turbulence 
kinetic  energy  for  this  obstacle  array.  In  addition  to  the  high-resolution  CFD 
simulation,  we  investigate  the  utility  of  neglecting  the  detailed  structure  of  the  solid 
boundaries  of  the  obstacles  in  the  array  and  considering  them  simply  as  an  aggregate 
effect  (viz.,  we  consider  the  aggregation  of  groups  of  buildings  in  the  array  and 
describe  them  simply  by  a  distributed  drag  force  on  the  airflow  entering  the  array). 


2.  The  wind  tunnel  experiment 

The  wind  tunnel  experiment  is  fully  described  in  Ref.  [8]  and  only  the  important 
details  of  the  experiment  will  be  presented  here.  The  experiments  were  conducted  in 
an  open-return  type  wind  tunnel  at  the  US  Environmental  Protection  Agency  s  Fluid 
Modeling  Facility.  This  wind  tunnel  has  a  working  test  section  of  length  18.3  m, 
width  3.7  m,  and  an  adjustable  roof  of  height  approximately  2.1m  to  eliminate 
streamwise  pressure  gradients  and  allow  for  a  non-accelerating  free  stream  flow. 

A  schematic  of  the  two-dimensional  building  array  used  for  these  wind  tunnel 
experiments  is  shown  in  Fig.  1.  This  array  consisted  of  seven  rectangular  blocks,  each 
block  with  equal  height,  H,  and  length,  L,  of  150  mm.  The  blocks  were  spaced  H 
apart  in  the  streamwise  direction.  The  building  array  was  immersed  in  a  simulated 
neutral  atmospheric  boundary  layer  which  was  created  in  the  wind  tunnel  using 
spires  and  floor  roughness  with  a  roughness  length  of  approximately  1  mm.  The 
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X/H 


I  ig  1  Computational  grid  superimposed  on  the  flow  domain  used  for  ihe  high-resolution  numerical 
simulation  of  the  two-dimensional  ( 2 1  >  >  array  of  buildings 


mean  streamwise  velocity  in  the  upstream  approach  flow  can  be  approximated  by  tlic 
following  power-law  form: 


where  u ti  3  m  s  1  is  the  reference  velocity  at  -  //  The  Reynolds  number*  based 

on  the  building  height  //  and  the  reference  velocity  hu,  is  30*000. 

The  streamwise,  cross-stream  (or,  span  wise),  and  vertical  velocity  components 
were  measured  with  a  pulscd-wire  anemometer  using  a  pulsing  rate  of  10  Hz  and  a 
sampling  time  of  120  s  at  each  measurement  location*  The  pulsed-wire  anemometer* 
unlike  a  standard  X -array  hot-wire  anemometer,  can  be  used  to  measure  velocities  in 
high  intensity  turbulence  with  reversing  flow's  such  as  those  which  exist  within  an 
obstacle  array.  Measurements  of  the  three  components  of  the  instantaneous  velocity 
were  made  at  1016  coordinate  positions.  These  measurements  consisted  of  81  vertical 
profiles  of  mean  velocity,  turbulence  velocity  variance,  and  turbulence  kinetic 
energy.  These  profiles  extended  from  3.5//  upstream  of  the  lirsi  building  to  7.3// 
downstream  of  the  last  building*  Each  vertical  profile  consisted  of  either  8  (over  the 
rooftops)  or  16  (over  the  urban  canyons,  and  either  upstream  or  downstream  of  the 
array)  measurement  positions  up  to  a  vertical  height  of  3//  above  the  ground 
surface.  All  streamwise  distances  arc  measured  relative  to  the  location  of  the 
upstream  wall  of  the  lirsi  building  at  xfii  0. 


3,  High-resolution  computational  Hu  id  dynamics 

3.  /.  Mathematical  model 

The  governing  equations  for  the  flow  of  an  incompressible  fluid  at  high  Reynolds 
number  \Re)  based  on  the  Reynolds-averaged  Navier  Stokes  approach  are: 
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5m,- 

dXj 


0, 


(2) 


5m,  Uj  _  1  dp  dt/jUj 
5  Ay  P  0X,  0Xy  ’ 


(3) 


where  the  Reynolds  averaging  of  a  quantity  is  denoted  by  drawing  a  bar  over  the 
quantity.  Here  u ,  and  mJ  are  the  mean  and  fluctuating  velocities  in  the  x, -direction, 
respectively;  p  is  the  density  of  the  fluid;  and,  p  is  the  dynamic  pressure.  In  Eqs.  (2) 
and  (3),  the  Einstein  summation  convention  is  assumed.  Subscripts  1-3  refer  to 
directions  along  x  (streamwise),  y  (spanwise),  and  z  (vertical),  respectively. 

Furthermore,  (u\ ,  U2,  M3)  =  ( u ,  v,  u  ).  _ 

The  presence  of  the  kinematic  Reynolds  stresses  m'mJ  in  the  mean  momentum 
equation  [Eq.  (3)]  implies  the  latter  is  not  closed.  Closure  requires  that  some 
approximation  be  made  in  prescribing  the  Reynolds  stresses  in  terms  of  the  mean 
flow  quantities.  A  common  closure  approximation  for  m'm',  first  suggested  by  Pope 
[9],  is 


m'm) 


5m,  duk  duj  duk\ 
dxk  dxj  +  dxk  dxj) 

"”V 

quadratic  term 


* 


/5 My-  5ma-\* 
\dxj  dxj) 


+  HOT, 


(4) 


where  5&  is  the  Kronecker  delta  symbol,  HOT  denotes  higher-order  terms,  and  l*' 
indicates  the  deviatoric  part;  for  example, 


/ dtij  duj \  5m,  duj  1  dum  5m„, 

\dxk  dxk  J  ~  dxk  dxk  3  5x„  5x„  v 


(5) 


Furthermore,  k  =  V'm-  is  the  turbulence  kinetic  energy  (TKE),  £  is  the  rate  of 
dissipation  of  turbulence  kinetic  energy,  and  vy  is  the  kinematic  eddy  (turbulent) 
viscosity.  In  the  conventional  linear  eddy-viscosity  models  (e.g.  Refs.  [10-12]),  only 
the  linear  term  is  retained.  For  Shih  et  al.’s  quadratic  non-linear  eddy-viscosity  model 
[13],  HOT  is  zero  and 


(Cr|,  Cr2,  Cxi) 


(■ 


13  -4  -2 _ \ 

1000  +  Sy  1000  +  S3’  1000  +  S3/  ’ 


(6) 


where 


S  = 


y/ISjjSjjk  c  _  1  (di!i  t  du/\ 
e  ’  *,j~  2\dxj  dx,)' 


(7) 


Here,  S  is  the  characteristic  strain  invariant  constructed  from  the  mean  rate-of-strain 
tensor  S,y. 
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In  the  framework  of  k  c  models,  the  eddy  viscosity  vy  is  given  by 
k2 

VT  =  c"  T’ 

where  k  and  e  are  determined  from  the  following  transport  equations: 
dbjk 


k_  0  (vT\dk' 
i  dxi[\'Tk)dxJ 

/ VT\  0£  ' 

\ajdxj_ 


diijS  _  6 
dxj  dxj 


+  Pk  ~  e, 


+  £  {CR\Pk  — 


with  Pk  the  production  of  TKE  being  defined  as 
*  - 


(8) 

(9) 

(10) 


(11) 


These  equations  contain  five  closure  constants;  namely,  C,„  ak,  at,  C£l,  and  Ct2. 
Different  variants  of  the  above  model  arise  from  the  different  approaches  for 
determining  the  values  for  these  coefficients.  Here,  four  different  high-/?e  k-e 
models  are  examined: 


1.  the  standard  (STD)  k-e  model  [10]; 

2.  the  Kato-Launder  (K-L)  k-e  model  [1 1]; 

3.  the  renormalization  group  (RNG)  k-e  model  [12]; 

4.  the  non-linear  (N-L)  eddy  viscosity  k-e  model  [13], 

The  first  three  models  are  based  on  the  linear  (Boussinesq)  stress-strain  relation. 
The  model  constants  associated  with  the  STD  and  K-L  turbulence  models  are 

Cn  =  0-09,  ak  =1,  ae=  1.3,  Cc,  =  1.44,  C,2  =  1.92.  (12) 

The  major  difference  between  the  STD  and  K-L  turbulence  models  lies  in  the 
modeling  of  the  TKE  production  (Pk)  term.  Instead  of  using  Eq.  (1 1)  as  in  the  case 
of  the  STD  model,  which  can  be  re-written  as 

Pk  =  Cf,eS 2,  (13) 

for  K-L  model, 

Pk  (14) 

where 


fl  =  yggg5fc  ifduL_duL\ 

e  2  \dxj  dxj ) 


(15) 


Here,  S2  is  the  characteristic  vorticity  invariant  that  is  constructed  from  the  mean 
rotation  tensor  Q, y.  The  K-L  turbulence  model  is  designed  specifically  to  suppress 
excessive  turbulence  energy  production  in  stagnation  regions.  In  particular,  for 
simple  shear  flows  SxQ  so  predictions  obtained  from  the  K-L  model  are  essentially 
identical  to  that  of  the  standard  model,  while  in  a  stagnation  region  £>=t:0  so  the 
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spurious  turbulence  energy  production  is  eliminated  here.  For  the  RNG  k-e  model, 


C„  =  0.085, 
CrA  =  1.44- 


<rk  =  0.72, 
5(1  -  S/4.38) 
1  +0.012S3 


ae.  —  0.72, 

,  Ce2  =  1.68. 


(16) 


Finally,  for  the  N-L  eddy  viscosity  model, 

„  0.667  ,  .  , 

C„  = - s - fffc  =  1,  ffe  =  1.3, 

'  1.25  +  5  +  0.9(2 

C,  =  1.44,  02  =  1.92. 


(17) 


3.2.  Numerical  framework 

Our  computational  domain  for  simulating  the  two-dimensional  building  array 
used  in  the  wind  tunnel  experiment  (cf.  Fig.  1)  extended  in  the  streamwise  direction 
from  x/H  -  -5  to  28,  with  the  upstream  wall  of  the  first  building  at  x/H  =  0.  The 
height  of  the  domain  was  9 H.  The  results  of  the  numerical  simulation  reported  below 
were  obtained  with  an  orthogonal  251  x  51  grid  shown  in  Fig.  1.  The  grid  lines  were 
concentrated  near  the  solid  surfaces  (ground,  rooftop,  and  building  walls)  and  the 
spacing  between  grid  lines  was  gently  stretched  at  increasing  distances  from  the  solid 
surfaces. 

Since  the  computational  domain  is  chosen  to  be  large  compared  with  the  two- 
dimensional  array  of  obstacles,  the  flow  at  its  boundary  is  not  affected  by  any  of  the 
obstacles.  Hence,  free  boundary  conditions  are  imposed  at  all  air-to-air  boundaries 
in  the  flow  domain.  At  the  inflow  boundary,  the  measured  distributions  of 
streamwise  mean  velocity,  w  =  u\,  the  vertical  mean  velocity,  tv  =  uj  =  0,  and  the 
turbulence  kinetic  energy,  k,  are  used.  This  information  is  obtained  from  the  vertical 
profiles  of  u  and  k  measured  furthest  upstream  from  the  upstream  face  of  the  first 
building  (viz.,  at  x/H  =  -3.5,  orx=  -500mm,  where  x  =  0  corresponds  to  the 
upstream  wall  of  the  first  building).  Above  z  =  3 H,  measurements  for  u  and  k  are 
not  available;  and,  for  the  inflow  boundary  conditions,  the  vertical  profile  of  u 
above  z  =  3H  was  assumed  to  adhere  to  a  power-law  profile  described  in  Eq.  (1), 
and  the  vertical  profile  of  k  was  assumed  to  be  constant  above  z  =  3 H 
(corresponding,  as  such,  to  a  deep  constant  stress  layer).  Unfortunately,  at  the 
inflow,  e  was  not  measured  and,  consequently,  was  estimated  simply  as  £  =  ki/2/l 
where  the  dissipation  length  scale  /  was  chosen  as  approximately  1.2 H.  The  length 
scale  /  is  unknown  at  the  inflow  boundary,  but  variation  of  its  value  from  0.1  H  to 
1.2 H  showed  that  it  had  only  a  very  small  impact  on  the  final  solution.  Far 
downstream,  at  the  outflow  boundary,  we  set  0w/0x  =  dk/dx  =  0e/0x  =  tv  =  0.  At 
the  upper  boundary,  we  set  du/dz  =  w  =  0,  dk/dz  =  ds/dz  =  0.  Finally,  at  all  solid 
boundaries  (ground,  obstacle  walls,  obstacle  roofs)  in  the  computational  domain, 
standard  wall  functions  are  used  for  mean  velocities  and  turbulence  quantities  [14]. 

The  numerical  method  employs  a  structured,  non-orthogonal  (allowing  boundary- 
fitted  coordinates),  fully  colocated  (viz.,  all  flow  variables  are  stored  at  the  same 
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location),  cell-centered,  finite  volume  approach  for  the  discretization  of  the 
computational  domain.  Within  this  non-orthogonal  system,  the  velocity  vector  is 
decomposed  into  its  Cartesian  components,  and  these  are  the  components  to  which 
the  momentum  equations  relate.  Advective  volume  face  fluxes  are  approximated 
using  a  third-order  upwind  discretization  scheme  [15],  and  physical  diffusive  cell 
fluxes  are  approximated  using  a  conventional  second-order  central  differencing 
scheme.  The  SIMPLE  algorithm  [16]  was  used  for  pressure  correction.  Here,  mass 
continuity  is  enforced  by  solving  a  pressure  correction  equation  which,  as  part  of  the 
iterative  sequence,  steers  the  pressure  toward  a  state  at  which  all  mass  residuals  in  the 
cells  are  negligibly  small.  In  conjunction  with  the  colocated  grid  used  here,  this 
method  is  known  to  provoke  chequerboard  oscillations  in  the  pressure  field, 
reflecting  a  state  of  pressure-velocity  decoupling.  To  avoid  this,  the  widely  used 
method  of  Rhie  and  Chow  [17]  is  adopted  to  non-linearly  interpolate  the  cell  face 
velocities  from  the  nodal  values  (at  the  center  of  the  cells).  The  interpolation  scheme 
essentially  introduces  a  fourth-order  pressure  diffusion  correction  that  smooths  out 
the  pressure  if  it  oscillates  rapidly.  In  the  numerical  scheme,  the  transport  equations 
for  the  mean  velocity,  turbulence  kinetic  energy,  and  viscous  dissipation  rate  and  the 
pressure-correction  equation  are  solved  sequentially  and  iterated  to  convergence, 
defined  by  reference  to  the  sum  of  the  absolute  cell  residuals  for  mass  and 
momentum  components.  A  tri-diagonal  matrix  algorithm  (TDMA)  applied 
iteratively,  in  a  line-by-line  fashion,  was  utilized  to  solve  the  system  of  algebraic 
equations  arising  from  the  discretization  of  the  transport  equations  and  pressure- 
correction  equation.  Iterative  refinement  of  all  fields  was  continued  until  the  sum  of 
the  absolute  cell  residuals  for  mass  and  for  momentum  was  satisfied  to  within  0.1% 
of  the  total  mass  and  momentum  fluxes,  respectively,  at  the  inflow  boundary. 
Further  details  of  the  numerical  procedure  can  be  found  in  Ref.  [18]. 

3.3.  Results  and  discussion 

Figs.  2  and  3  show  the  mean  flow  streamlines  superimposed  on  the  turbulence 
kinetic  energy  field  in  the  x-z  plane  around  and  upstream  of  the  first  three 
buildings  and  around  and  downstream  of  the  last  building  in  the  array,  respectively, 
for  the  standard,  K-L,  RNG,  and  non-linear  k-e  turbulence  models.  Qualitatively, 
a  number  of  gross  features  in  the  flow  are  correctly  predicted  by  these  models.  The 
separation  zone  where  the  flow  detaches  from  the  windward  rooftop  edge  of  the  first 
building  and  reattaches  to  the  roof  near  the  leeward  rooftop  edge,  and  the  absence  of 
such  a  separation  zone  on  the  subsequent  (downwind)  rooftops  is  well  predicted,  as 
is  the  single  recirculating  vortex  that  forms  in  the  street  canyon  between  the 
buildings.  The  stagnation  point  on  the  windward  face  of  the  first  building  which 
occurs  at  about  one-half  the  building  height  is  correctly  predicted  by  the  various 
models.  The  slow-moving  and  relatively  large  recirculation  zone  formed  behind  the 
leeward  wall  of  the  last  building  has  been  predicted,  although  two  of  the  turbulence 
models  (namely,  the  RNG  and  non-linear  models)  appear  to  overpredict  and  the 
other  two  models  (namely,  the  standard  and  K-L  models)  appear  to  underpredict 
(albeit  slightly)  the  experimental  reattachment  length  (%3.8 H)  of  this  separation 
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zone.  All  the  turbulence  models,  with  perhaps  the  exception  of  the  K-L  model 
(which  was  designed  to  solve  the  stagnation  point  anomaly  by  reducing  the 
production  of  TKE  in  regions  of  high  streamwise  strain),  exhibit  a  strong  gradient  of 
TKE  centered  near  the  top  and  slightly  upstream  of  the  windward  edge  of  the 
rooftop  of  the  first  building  with  the  TKE  “exported”  to  locations  downstream  of 
this  position  in  a  broad  TKE  “plume”. 

3.3.1.  Mean  streamwise  velocity 

Fig.  4  compares  model  predicted  vertical  profiles  of  the  mean  streamwise  velocity 
u  at  four  fixed  x-locations  covering  the  range  from  within  the  impingement  zone 
directly  upstream  of  the  first  building,  over  the  rooftop  of  the  first  building,  into  the 
first  street  canyon,  and  over  the  rooftop  of  the  second  building.  At  x  =  —75  mm,  the 
wind  speed  reduction  below  and  wind  speed  increase  above  the  level  z/Hx 4/3 
relative  to  the  upstream  profile  is  reproduced  well.  However,  the  models  do  not 
predict  the  reversal  in  the  mean  flow  close  to  the  ground  at  this  position. 
Furthermore,  close  to  the  upstream  wall  of  the  first  building  at  x=  -15  mm  (not 
shown),  the  significant  wind  speed  reduction  below  the  building  height  H  owing  to 
the  adverse  pressure  rise  (horizontal  pressure  gradient)  upstream  of  the  solid  obstacle 
is  modeled  well,  as  is  the  moderate  speed-up  above  the  building  rooftop.  At 
x  =  75  mm,  the  speed-up  over  the  first  rooftop  compared  with  the  reference  upwind 
profile  and  the  presence  of  a  thin  separation  zone  near  the  roof  surface  with  flow 
reversal  are  well  matched  between  the  predictions  and  the  measurements.  The  results 
for  .v  =  225  mm  show  reasonable  agreement  between  the  measurements  and  the 
model  predictions  of  u  in  the  first  urban  canyon.  An  important  feature  captured  by 
the  model  predictions  here  is  the  very  strong  shear  layer  that  forms  immediately 
downstream  of  the  rear  face  of  the  first  building,  whose  signature  is  revealed  by  the 
inflection  point  in  u(z)  at  or  near  the  building  height  H.  Note  that  the  high  values  of 
mean  shear  9 u/dz  just  above  building  height  over  the  first  urban  canyon  are 
reproduced  well  by  the  model.  The  speed-up  above  about  z/Hx 4/3  is  predicted 
correctly  over  the  first  urban  canyon,  although  the  magnitude  of  the  reverse  flow  in 
the  canyon  below  about  z///~  2/3  is  underestimated  slightly  by  the  standard  k— e 
model.  The  reverse  velocity  on  the  underside  of  the  standing  vortex  within  the 
canyon  is  comparable  to  the  forward  velocity  at  the  canyon  top,  and  this  feature 
appears  to  be  reproduced  by  the  model  predictions.  Finally,  the  computed  and 
measured  mean  streamwise  velocity  u  over  the  rooftop  of  the  second  building  is  well 
predicted.  The  models  predicts  that  no  separation  occurs  over  the  second  rooftop  (in 
contrast  to  the  behavior  over  the  first  rooftop)  which  accords  well  with  the 
observations. 

Vertical  profiles  of  w  at  four  x-locations  within  the  second  and  third  urban 
canyons  and  over  the  rooftops  of  the  third  and  fourth  buildings  are  displayed  in 
Fig.  5.  Note  that  the  predicted  magnitudes  of  the  wind  shear  over  the  second  and 
third  street  canyons  are  larger  than  what  is  observed  in  the  wind  tunnel 
measurements.  The  primary  reason  for  this  discrepancy  appears  to  be  the  reduced 
vertical  turbulent  diffusion  of  shear  stress  in  the  high  shear  layer  bordering  the  urban 
canyon  top— a  defect  manifesting  itself  by  the  greater  level  of  shear  strain  in  this  area 
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than  observed.  The  reduced  vertical  turbulent  diffusion  results  from  the  under¬ 
estimation  of  the  turbulence  kinetic  energy,  which  implies  an  underestimation  in  the 
turbulent  viscosity.  Because  of  this,  the  high  shear  layer  generated  near  the  urban 
canyon  top  dissipates  more  slowly  in  the  streamwise  direction  along  the  canyon  and 
the  succeeding  rooftops  than  in  the  wind  tunnel  observations.  In  consequence,  the 
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Fig.  5.  Vertical  profiles  of  the  mean  streamwise  velocity  u  at  four  .v-locations  (.v  =  525,  675,  825  and 
975  mm),  obtained  from  a  high-resolution  numerical  simulation  using  four  different  turbulence  models, 
are  compared  with  time-averaged  wind  tunnel  measurements  at  the  same  locations. 


prediction  of  the  speed-up  above  and  wind  reduction  below  the  level  z/H^ 4/3  is 
slightly  over-predicted  and  under-predicted,  respectively. 

Vertical  profiles  of  both  the  computed  and  measured  mean  streamwise  velocity  u 
within  the  third  urban  canyon  and  the  fourth  rooftop  are  almost  identical  to  the 
profiles  obtained,  respectively,  for  the  second  urban  canyon  and  the  third  rooftop.  In 
consequence,  u  appears  to  have  reached  streamwise  equilibrium  (i.e.,  the  mean 
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Fig.  7.  Vertical  profiles  of  the  turbulence  kinetic  energy  k  at  four  v-locations  (.x:  =  -75,  75,  225  and 
375  mm),  obtained  from  a  high-resolution  numerical  simulation  using  four  different  turbulence  models, 
are  compared  with  time-averaged  wind  tunnel  measurements  at  the  same  locations. 


zone  would  correspond  to  a  gain  in  -mV  owing  to  the  positive  (unstable)  flow 
curvature  (/?>0)  in  this  region  (noting  that  2m'2S>h’/2).  Hence,  a  reduction  in 
normal-stress  anisotropy  would  lead  to  an  underestimation  in  -iTw7,  which  in  turn 
would  lead  to  a  reduced  shear  stress  production,  -uV8fl/3z,  of  the  TKE.  Secondly, 
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an  underestimation  of  the  normal  stress  anisotropy  will  lead  to  an  under-prediction 
of  the  normal  stress  production  of  k ,  ~{ua  -  h’/2)Sw/9x,  in  the  impingement  zone  in 
front  of  the  first  building.  Note  that  in  this  region  du/dx<0  owing  to  the 
deceleration  of  the  flow  in  front  of  the  windward  face  of  the  first  building  (induced 
by  the  adverse  horizontal  pressure  gradient  here),  so  the  term  (in  streamline 
coordinate  system)  —  (a'2  —  \vf2)du/dx,  corresponds  to  a  gain  in  TKE. 

If  we  approximate  ~(u'2  —  wf2)du/dx  by  4 vr(du/dx)2  following  the  Boussinesq 
linear  stress-strain  relation,  then  an  excessive  turbulence  kinetic  energy  is  expected  to 
be  generated  along  the  stagnation  line  due  to  sufficiently  high  normal  straining — a 
defect  commonly  referred  to  as  the  ‘stagnation  point  anomaly’.  This  results  in  the 
standard  k-e  model  overestimating  the  TKE  levels  in  the  stagnation  zone  of  a  bluff 
body  (see,  for  example,  [14]).  However,  this  seems  to  contradict  the  observed 
underestimation  of  the  predicted  levels  of  TKE  in  the  impingement  zone  by  the 
standard  k-e  model  for  the  current  situation.  A  possible  explanation  for  this 
contradiction  is  that  the  standard  k-e  model  does  not  account  properly  for  the 
subtle  effects  of  streamline  curvature  on  the  turbulence.  In  the  present  case,  the 
impingement  process  occurring  in  front  of  the  first  building  wall  provokes  a  severe 
concave  (unstable)  flow  curvature  in  the  region  of  the  flow  immediately  upstream  of 
the  first  building  rooftop.  It  is  expected  that  the  turbulence  kinetic  energy  in  this 
region  will  increase  significantly  due  to  the  cumulative  effect  of  concave  streamline 
curvature  as  fluid  elements  are  forced  up  and  over  the  first  building.  The  curvature 
term  (in  streamline  coordinates)  -u'w'(u/R)  is  a  source  of  production  of  k  in  this 
region  of  positive  (unstable)  R ,  and  this  source  of  TKE  is  not  properly  modeled  in 
the  standard  k-e  model. 

Alternatively,  the  observed  large  peak  in  the  TKE  just  above  the  first  building 
rooftop  can  be  also  interpreted  as  arising  from  an  oscillation  of  the  high  shear  layer 
by  the  larger-scale  upstream  turbulence  which  would  enhance  u'2  (and,  hence  k). 
Certainly,  this  “flapping”  of  the  shear  layer  at  the  building  height  over  a  street 
canyon  has  been  reported  in  Ref.  [19].  The  fact  that  such  a  large  peak  is  missing  in 
the  predictions  may  suggest  that  this  process  is  not  well  represented  in  the  k—e 
turbulence  model.  Indeed,  steady  RANS  predictions  cannot,  of  course,  account  for 
these  large-scale  unsteady  motions  (i.e.,  shear  layer  “flapping”).  Since  TKE  is 
exported  downstream  by  local  advection  and  turbulent  transport,  an  under¬ 
estimation  of  the  TKE  level  in  the  region  just  upstream  of  the  first  building  will 
lead  generally  to  an  underestimation  of  the  TKE  level  at  all  positions  downstream  of 
the  location  of  strong  turbulence  generation,  although  this  defect  is  expected  to 
diminish  with  increasing  downstream  distance  from  the  region  of  maximum  TKE 
production. 

The  position  of  the  prominent  “nose”  in  the  k  profiles  lying  just  above  the  street 
canyon  top  at  x  =  225  mm  is  largely  reproduced  by  the  numerical  simulation 
(especially  by  the  standard  model),  albeit  the  peak  value  of  k  is  underestimated.  This 
defect  is  a  consequence  of  the  under-prediction  of  the  TKE  level  upstream  of  this 
region  owing  to  the  inability  of  the  present  models  to  account  properly  either  for  the 
effects  of  secondary  strain  on  the  turbulence  and/or  for  the  effects  of  the  large-scale 
flapping  of  the  strong  shear  layer  at  the  canopy  top.  Moving  downstream  from  this 


134 


F.S  Lien  et  al  f  J.  Wind  Eng.  Ind  Aerodyn.  92  (2004)  117-158 


position,  the  evolution  of  the  TKE  profiles  is  dominated  by  the  vertical  spreading  of 
this  shear  or  mixing  layer  (cf.  Fig.  8).  At  the  center  of  this  layer  lying  just  above  the 
building  height  H ,  the  peak  value  of  TKE  attenuates  downstream.  Turbulence 
energy  is  exported  to  regions  below  and  above  the  peak  by  pressure  or  turbulent 
transport  from  the  shear  layer  core,  causing  the  TKE  to  increase  in  these  regions. 


Fig.  8.  Vertical  profiles  of  the  turbulence  kinetic  energy  k  at  four  ^locations  (.v  =  525,  675,  825  and 
975  mm),  obtained  from  a  high-resolution  numerical  simulation  using  four  different  turbulence  models, 
are  compared  with  time-averaged  wind  tunnel  measurements  at  the  same  locations. 
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These  characteristics  of  the  flow  are  largely  reproduced  by  the  numerical  simulations 
(especially  by  the  standard  and  N-L  turbulence  models). 

The  downstream  advection  and  vertical  diffusion  of  TKE  results  in  a  larger  level 
of  turbulence  energy  being  found  on  the  downstream  (leeward)  side  of  a  street 
canyon  than  the  upstream  (windward)  side,  a  feature  that  is  also  present  in  the  model 
predictions  (not  shown).  Within  the  street  canyon,  minimum  levels  of  TKE  are 
found  on  the  upstream  side  of  the  canyon  (protected  region)  where  the  TKE  is 
essentially  independent  of  height.  These  features  are  reproduced  in  the  simulations, 
although  the  near  constant  TKE  level  on  the  upstream  side  of  the  canyon  is 
underestimated  by  the  models.  In  addition,  the  TKE  levels  in  street  canyons  further 
downstream  are  successively  smaller  owing  to  the  fact  that  the  peak  TKE  in  the 
mixing  layer  above  each  canyon  attenuates  downstream  (and,  hence,  there  is  less 
TKE  available  to  transport  into  each  successive  street  canyon  downstream).  These 
qualitative  features  of  the  flow  are  largely  reproduced  by  the  numerical  simulations. 
The  quantitative  agreement  between  the  computed  and  measured  TKE  profiles 
generally  improves  with  increasing  downstream  distance.  Finally,  it  appears  that  the 
rate  of  vertical  spreading  of  the  shear  or  mixing  layer  in  the  numerical  simulation  is 
smaller  than  in  the  measurements.  This  discrepancy  may  be  the  result  of  an 
underestimation  of  the  turbulent  viscosity  vr  owing  to  the  smaller  predicted  TKE 
levels. 

Fig.  9  compares  computed  and  measured  TKE  profiles  near  the  exit  region  of  the 
array.  Overall,  the  TKE  profiles  in  the  exit  region  appear  to  be  very  similar  in  shape, 
and  all  the  qualitative  features  in  these  profiles  appear  to  be  correctly  predicted;  the 
quantitative  agreement  is  also  fairly  good,  except  for  the  RNG  model.  For  the 
profiles  in  the  range  from  2025  to  2775  mm  (only  profiles  at  x  =  2175  and  2475  mm 
are  shown  here  and  both  locations  are  within  the  recirculation  zone),  the  position  of 
the  maximum  TKE  is  fairly  well  predicted  by  all  models.  This  position  corresponds 
to  the  center  of  the  curved  shear  layer  that  bounds  or  envelops  the  recirculation  zone 
behind  the  last  building.  Finally,  the  rate  of  recovery  of  the  TKE  in  the  far  wake 
region  to  the  upstream  reference  condition  is  reproduced  correctly.  Even  at 
*  =  3075  mm  or,  equivalently,  x/H  =  1.5  downstream  of  the  last  building 
(measurement  furthest  downstream),  the  TKE  profile  has  the  appearance  of  a 
detached  mixing  layer  despite  the  fact  that  the  mean  flow  has  re-attached  at  about 
x/H  =  3.8  downstream  of  the  last  building  (not  shown).  This  profile  form  is  a 
residue  of  the  strong  shear  or  mixing  layer  that  detaches  from  the  last  building 
rooftop  and  spreads  outward  by  pressure  or  turbulent  diffusion  until  this  vertically 
spreading  mixing  layer  makes  contact  with  the  ground  surface.  This  is  an  example  of 
disequilibrium  that  results  in  the  TKE  profile  at  even  the  furthest  downstream 
measurement  position  not  having  recovered  to  the  reference  far  upstream  state. 


3.3.3.  Summary  of  model  performance 

Based  on  the  preceding  results,  it  is  found  that 

(1)  the  non-linear  k~e  model  is  capable  of  overcoming  the  deficiency  of  the 
stagnation  point  anomaly  associated  with  the  standard  k~e  model  and. 
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Fig.  9.  Vertical  profiles  of  the  turbulence  kinetic  energy  k  at  four  ^-locations  (x  =  1725,  1875,  2175  and 
2475  mm),  obtained  from  a  high -resolution  numerical  simulation  using  four  different  turbulence  models, 
are  compared  with  time-averaged  wind  tunnel  measurements  at  the  same  locations. 


therefore,  gives  a  better  prediction  of  the  mean  streamwise  velocity  profile  in  the 
street  canyon  at  x  =  225  m  (cf.  Fig.  4); 

(2)  the  non-linear  k~e  model  provides  better  predictions  of  the  mean  streamwise 
velocity  profiles  at  x  =  2175  and  2475  mm  in  the  wake  region  behind  the  last 
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building  (cf.  Fig.  6),  especially  above  the  canopy  height,  than  the  standard  k-e 
and  K-L  turbulence  models,  although  it  overestimates  the  reattachment  length; 
(3)  the  non-linear  k-e  model  provides  better  predictions  of  the  turbulence-energy 
profiles  within  the  street  canyon,  above  the  rooftops  of  buildings,  and  in  the 
wake  region  downstream  of  the  last  building  compared  to  the  RNG  k-e  model. 

In  addition,  the  non-linear  k-e  model  is  algorithmically  simpler  and  numerically 
more  robust  and  efficient  than  second-moment  closure  models.  However,  in  many 
cases,  the  non-linear  k-e  model  has  a  similar  predictive  performance  compared  with 
the  (necessarily  more  complex)  second-moment  closure  models  [20,21].  Therefore,  we 
recommend  that  the  non-linear  k-e  model  be  used  for  the  prediction  of  flows  over  a 
cluster  of  buildings. 


4.  Distributed  drag  force  approach 

4.1.  Mathematical  model 

In  the  previous  section,  all  the  buildings  in  the  array  were  resolved  explicitly  in  the 
sense  that  boundary  conditions  were  imposed  at  all  the  walls  and  roofs  of  the 
buildings.  In  this  section,  we  will  investigate  the  utility  of  representing  individual 
buildings  and/or  groups  of  buildings  in  the  obstacle  array  in  terms  of  a  distributed 
drag  force.  The  process  will  involve  the  decomposition  of  the  obstacle  array  into 
various  contiguous  regions  that  contain  one  or  more  obstacles  with  the  aggregate 
effect  of  the  obstacle(s)  on  flow  in  these  regions  represented  by  a  drag  coefficient. 
This  approach  will  obviate  the  need  to  explicitly  impose  boundary  conditions  on  the 
surfaces  of  all  buildings  (and  other  obstacles)  in  the  array.  A  similar  approach  has 
been  used  by  Belcher,  Jerram,  and  Hunt  [22]  to  study  the  adjustment  of  the  mean 
velocity  to  a  canopy  of  roughness  elements  using  a  linearized  flow  model  (obtained 
by  determining  analytically  small  perturbations  to  the  undisturbed  upstream 
logarithmic  mean  velocity  profile  induced  by  the  drag  due  to  an  obstacle  canopy). 
However,  this  study  investigated  only  the  adjustment  of  the  spatially  averaged  time- 
mean  flow  over  a  roughness  element  canopy,  and  did  not  focus  on  any  of  the 
turbulence  statistics  (e.g.,  turbulence  kinetic  energy). 

The  inclusion  of  additional  source/sink  terms  in  the  mean  momentum  equations 
and  the  supporting  k -  and  e-equations  to  represent  the  aggregate  effect  of  the 
obstacles  on  the  flow  can  be  obtained  if  we  consider  first  the  addition  of  an 
instantaneous  (i.e.,  before  the  application  of  Reynolds  averaging)  drag  force  term 

/d,  =  CDA(UjUj)'/2uh  (18) 

which  removes  momentum  in  proportion  to  the  instantaneous  velocity,  in  the  un¬ 
averaged  Navier-Stokes  equation.  In  Eq.  (18),  «,  is  the  instantaneous  spatially 
averaged  velocity,  CD  is  the  drag  coefficient,  and  A  is  the  frontal  area  density  (i.e.,  the 
frontal  area  of  the  obstacles  per  unit  volume  of  space).  Now,  using  the  Reynolds 
decomposition  to  decompose  u ,  into  its  mean  u,  and  fluctuating  ili  components  (viz.. 
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u,  =  tij  +  wj),  the  instantaneous  drag  force  can  be  expanded  and  then  approximated 
up  to  third-order  as  follows: 


/d  /  =  CDA[(Uj  +  Uj)(tij  +  u'j)]'/2(u,  +  u'i) 


CdA 


=  CdA 


i  n  ui^i 
(UjUj)1  + 


( akUk)l/ 2 


(«< + 


I  n  ,  UjUiU:  UjU’.u': 

(UjUj)  /  («,  +  M()  +  _  |/2  +  _  j  /, 

' - - - '  ( Ukuk )  7  ( UkUkY' 


III 


(19) 


In  the  second  line  of  Eq.  (19),  it  is  assumed  that  UjUjPu’jUj  and  the  binomial  theorem 
is  used  to  approximate  the  resulting  square  root  term  (UjUj  +  2uiu,j)l/2. 

The  three  terms  (namely,  I,  II,  and  III)  on  the  right-hand  side  of  Eq.  (19) 
correspond  to  various  levels  of  approximation  that  have  been  applied  to  the 
expansion  of  the  instantaneous  drag  force  term  in  order  to  facilitate  the  subsequent 
Reynolds  averaging  process  required  to  include  this  term  in  the  mean  momentum 
equations  (and,  also  in  the  supporting  transport  equations  for  k  and  e).  For  example, 
Getachew  et  al.  [23]  retain  terms  I,  II,  and  III  in  their  third-order  approximation  of 
the  instantaneous  drag  force  term.  However,  the  treatment  by  Antohe  and  Lage  [24] 
and  Ayotte  et  al.  [25]  propose  the  retention  of  terms  I  and  II  only  to  give  a  second- 
order  approximation  for  the  instantaneous  drag  force  term.  Finally,  Lee  and  Howell 
[28]  use  only  term  I  in  their  first-order  approximation  of  the  instantaneous  drag  force 
term.  We  have  used  the  drag  parameter  CqA  in  our  current  parameterization  of  the 
instantaneous  drag  force  term  [cf.  Eq.  (18)],  but  in  the  porous  media  literature  (e.g.. 
Refs.  [23,24])  this  drag  parameter  is  normally  replaced  by  F/k^2,  where  F  is  the 
Forchheimer  constant  (which  has  values  approximately  in  the  range  from  0.075  to 
0.5  according  to  measurements)  and  k  is  the  permeability.  It  is  noted  that  ki/2 
corresponds  to  a  length  scale  that  is  characteristic  of  the  effective  pore  size  of  the 
porous  medium.  Note  that  the  dimension  of  k,/2  is  the  same  as  that  of  A~l  (i.e.,  both 
are  expressed  in  units  of  length).  Hence,  it  is  possible  to  make  the  following 
identification:  CdA<->F/k'/2. 

If  we  ensemble-average  (or,  time-average)  the  instantaneous  drag  force  term  of 
Eq.  (19)  with  the  retention  of  either  terms  I  and/or  II  and  ignore  the  ‘dispersive 
stress’  term1  (see,  for  example,  [25])  for  simplicity  (Ockham’s  razor),  we  find  that  the 
drag  force  contribution  to  the  mean  momentum  equation  (which  is  the  prognostic 


1  The  'dispersive  stress'  corresponds  to  a  ‘dispersive’  momentum  flux  term  that  physically  arises  from  the 
correlations  in  the  point-to-point  variations  in  the  time-averaged  (mean)  velocity  field.  Mathematically, 
this  term  arises  from  spatially  averaging  the  non-linear  convective  term  in  the  mean  momentum  (RANS) 
equation.  The  sparse  experimental  evidence  to  date  seems  to  suggest  that  the  dispersive  stress  term  is 
negligible  in  comparison  to  the  spatially  averaged  Reynolds  stress  [26]  in  a  uniform  plant  canopy,  but 
whether  this  is  true  for  an  urban-type  roughness  array  remains  to  be  demonstrated  [27],  In  a  later  section, 
we  will  use  the  high-resolution  CFD  results  described  in  Section  3  to  diagnose  the  dispersive  stresses  for 
2D  building  array. 
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equation  allowing  the  prediction  of  the  spatially  averaged  time-mean  velocity)  has 
the  following  form: 

8 UjUj  _  1  dp  8m/“/  ?  ^  (20) 

8xj  pdxj  dxj  D’’ 


where 

/d,  =  CDA(UjUj)'/2Ui  (2I) 

is  the  mean  (ensemble-averaged)  drag  force  exerted  on  a  unit  mass  of  fluid  (or, 
specific  drag  force)  in  the  x,- -direction  induced  by  either  a  single  canopy  element  or  a 
group  of  canopy  elements.  In  this  section,  unless  otherwise  noted,  time  and  spatially 
averaged  velocity  components  (subscripted)  will  be  represented  with  overbars. 
Consequently,  in  Eq.  (20),  it,  is  interpreted  as  the  spatially  averaged  time-mean 
velocity  component  in  the  x, -direction  within  the  context  of  the  distributed  drag 
force  approach.  Furthermore,  in  Eq.  (20),  it  is  noted  that  /  />,  is  defined  to  be  positive 
in  the  direction  opposite  to  n,  (viz.,  the  flow  is  decelerated  if  f  d,  is  positive;  or, 
equivalently,  the  drag  force  always  acts  against  the  flow  direction).  The  term  f  D, 
represents  the  drag  of  the  obstacles  on  the  mean  momentum.  Finally,  it  is  noted  that 
a  third-order  approximation  to  the  instantaneous  drag  force  term  results,  on 
ensemble-averaging,  in  an  additional  term  in  the  mean  drag  force  involving  the 
Reynolds  stresses  [23]. 

So  far,  we  have  considered  only  the  addition  of  a  mean  drag  force  term  /  p,  to  the 
mean  momentum  equation.  The  effect  of  form  and  viscous  drag  in  the  urban 
canopy  on  the  k-  and  e-equations  will  also  require  the  inclusion  of  additional  source 
and  sink  terms  in  these  equations.  It  is  remarked  that  in  many  simulations  of  canopy 
flow  (or,  of  flow  in  porous  media),  no  additional  source  and  sink  terms  are  included 
in  the  k-  and  e-equations  (e.g.,  Refs  [29,30]),  but  as  noted  by  Wilson  [31]  with 
reference  to  a  second-order  closure  model  the  “quiet  zone’’  corresponding  to  the 
region  of  reduced  turbulence  in  the  lee  of  a  porous  barrier  cannot  be  properly 
simulated  without  the  inclusion  of  an  appropriate  sink  for  n'2  at  the  barrier.  In  view 
of  this,  we  will  adopt  the  idea  of  Lee  and  Howell  [28]  who  proposed  a  simple 
treatment  of  source  and  sink  terms  in  the  k— e  turbulence  model  for  the  case  of  a 
highly  porous  layer  exposed  to  a  turbulent  external  flow  field. 

In  the  present  study,  we  only  consider  Lee  and  Howell’s  approach  for  simplicity. 
The  other  more  complicated  approaches  (e.g..  Refs  [23,24])  lead  to  the  addition  of 
many  more  source/sink  terms  in  the  mean  momentum  and  k-  and  e-equations,  but 
their  accuracy  and  generality  have  not  been  proven;  we  prefer  to  retain  the  simplest 
model  until  the  need  for  a  better  representation  is  demonstrated.  In  consequence, 
using  the  Lee  and  Howell  [28]  approximation  for  the  instantaneous  drag  term  [e.g., 
term  I  in  Eq.  (19)],  the  evolution  equation  for  becomes 


0/ 


-  CoAiUjUj)'^. 


(22) 


Multiplying  Eq.  (22)  by  w',  applying  the  summation  in  i,  averaging  in  time 
(Reynolds  averaging),  and  using  the  definition  k  =  u '.«{/2  yields  the  evolution 
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equation  for  the  (spatially  averaged)  TKE: 

dk  *  m 

a7=-  -2CDA(UjUj)'/2k  . 

extra  sink  term  due  to  drag 

In  view  of  this,  the  evolution  equation  for  k  can  be  written  explicitly  as 
diijk  6  f  vy\  0&1 

Qxj  ~  Ws^J  +  {Pk  ~ £)  “  2c»A(m)  k- 


(23) 


(24) 


The  extra  sink  term  in  the  transport  equation  for  k  represents  the  loss  from  k  to 
smaller  scales  and  physically  corresponds  to  the  loss  of  energy  from  the  turbulent 
eddies  in  the  canopy  as  these  eddies  do  work  against  the  form  and  viscous  drag  of  the 
canopy  elements.  The  modeled  transport  equation  for  e  cannot  be  derived 
systematically.  However,  a  dimensionally  consistent  analogy  to  the  above  k- 
equation  can  be  written  down  as  a  modeled  evolution  equation  for  e\ 


4.2.  Determination  of  the  drag  coefficient 

The  approach  for  parameterizing  the  mean  drag  force  in  terms  of  CD,  described  in 
Section  4.1,  has  been  applied  previously  to  the  numerical  modeling  of  flows  in 
both  vegetative  (e.g.,  Refs.  [25,32-36];  among  others)  and  urban  (e.g.,  Refs.  [37-39]) 
canopies.  A  major  problem  with  this  approach  resides  in  the  assignment  of  an 
appropriate  value  for  CD  for  a  canopy  consisting  of  an  arbitrary  geometrical 
arrangement  of  elements.  This  is  a  non-trivial  and  unsolved  problem  for 
there  is  currently  no  general  method  that  can  be  used  to  predict  the  in  situ  drag 
coefficient  of  an  element  within  a  canopy  from  a  knowledge  of  the  drag  coefficient  of 

the  element  in  isolation  and  the  geometrical  arrangement  of  the  elements  within  the 
canopy. 

In  view  of  this  difficulty,  the  drag  coefficient  CD  has  usually  been  treated  as  an 
adjustable  model  parameter  whose  value  is  chosen  to  optimize  the  fit  between  the 
predicted  and  measured  mean  velocity  profiles  (e.g.,  Refs.  [25,33]).  Alternatively,  for 
horizontally  homogeneous  canopies  with  no  mean  horizontal  pressure  gradient, 
eithei  a  local  or  bulk  drag  coefficient  for  the  canopy  can  be  diagnosed  from  the 
measured  (observed)  profiles  of  shear  stress  and  mean  velocity  [40], 

For  a  horizontally  homogeneous  canopy  with  no  appreciable  streamwise  mean 
pressure  gradient  (i.e.,  for  fully  developed  flow),  the  mean  momentum  equation 
within  the  canopy  simply  involves  a  balance  between  the  stress  gradient  force  (i  e 
the  vertical  gradient  of  the  shear  stress)  and  the  drag  force  (imposed  on  the  flow  by 
the  canopy).  Hence,  the  measured  fi(z)  and  kinematic  shear  stress  f(z)  =  -w'tv'(z)  can 
be  used  to  diagnose  the  local  drag  coefficient  parameter  as 

coto-fg)*. 
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If  vertical  profiles  of  the  shear  stress  are  not  available,  a  bulk  drag  coefficient 
parameter  CDA  (which  is  assumed  to  be  constant  with  height)  can  be  determined  as 
follows: 

f (H)  =  CdA  [  u2  d z,  (27) 

Jo 

where  x(H)  is  the  kinematic  shear  stress  at  the  canopy  top  (i.e.,  at  z  =  H  where  H  is 
the  canopy  height).  Eq.  (27)  is  obtained  by  integrating  Eq.  (26)  with  respect  to  z  from 
ground  level  to  the  canopy  top  and  assuming  that  the  shear  stress  vanishes  at  the 
ground  (viz.,  f(z  =  0)  =  0;  or,  equivalently,  that  the  momentum  absorbed  by  the 
ground  is  negligible  compared  to  that  absorbed  by  the  canopy).  The  value  of  CDA 
determined  from  either  Eqs.  (26)  or  (27)  can  then  be  used  in  a  numerical  model  to 
predict  the  canopy  flow.  However,  given  the  assumptions  used  to  derive  either 
Eqs.  (26)  or  (27),  this  approach  is  not  applicable  for  developing  and/or  horizontally 
inhomogeneous  canopy  flows.  We  require  an  alternative  methodology  for  the 
determination  of  Co  in  this  case. 

The  drag  force  FDl  exerted  on  a  canopy  element  or  group  of  canopy  elements  in  a 
control  (averaging)  volume  V  in  the  x, -direction  is  the  surface  integral  of  the 
pressure  (form)  and  viscous  (frictional)  forces  over  all  canopy  element  surfaces 
contained  in  V  [25];  namely. 


where  n,  is  the  unit  normal  vector  in  the  x, -direction  pointing  away  from  the  surface 
of  the  canopy  element  into  the  fluid,  and  x y,-  is  the  dynamic  shear  stress  in  the  x,- 
direction  exerted  on  a  surface  (plane)  of  a  canopy  element  perpendicular  to  the  xy- 
direction.  Here,  the  shear  stress  xy,  on  the  surface  of  a  canopy  element  will  be 
obtained  explicitly  from  a  high-resolution  CFD  simulation  using  the  ‘wall-function’ 
approach  (e.g.,  Ref.  [14]).  Once  the  drag  force  fD,  has  been  determined,  a  drag 
coefficient  Co  can  be  diagnosed  using  this  information. 

In  this  paper,  we  will  explicitly  diagnose  a  drag  coefficient  CD  for  the  case  of  a 
developing  flow  over  an  array  of  2D  obstacles  using  high-resolution  CFD.  More 
specifically,  the  test  problem  here  is  the  same  as  in  Section  3.3  namely,  a  developing 
flow  over  an  array  of  seven  2D  buildings.  In  this  case,  Co  must  necessarily  account 
for  the  effects  of  the  drag  on  the  “urban”  canopy  turbulence  and  the  ‘sheltering’ 
within  the  canopy.  In  consequence,  Cd  does  not  necessarily  have  to  be  constant  for 
each  building  in  this  array.  Specifically,  it  is  anticipated  that  the  form  (pressure)  drag 
will  be  large  for  the  first  and/or  second  buildings  of  the  array  in  the  upstream 
direction  where  the  incident  flow  first  impinges  on  the  array.  To  proceed  with  the 
analysis,  we  decompose  the  2D  building  array  into  six  drag  units  as  illustrated  in 
Fig.  10(a).  Each  drag  unit  consists  of  one  building  and  the  associated  downstream 
street  canyon.  The  seventh  building  in  the  array  will  be  explicitly  resolved  (i.e., 
appropriate  boundary  conditions  will  be  imposed  along  the  walls  and  rooftop  of  this 
last  building). 


6  Drag  Units 


Pig.  10  Decomposition  of  the  array  into  fa)  six  drag  units  with  the  last  obstacle  of  the  array  explicitly 
resolved;  (b)  one  drag  unit  with  the  first  and  last  obstacles  of  the  array  explicitly  resolved;  and.  (c) 
expanded  view  of  one  drag  unit  from  (a)  showing  the  "Forces’"  acting  on  the  surface  of  the  drag  unit. 


Now,  consider  the  expanded  view  of  one  of  these  drag  units  shown  in  Fig.  10(c). 
The  control  volume  occupied  by  the  drag  unit  is  assumed  to  have  unit  length  in  the 
spanwise  (or.  y~)  direction.  The  balance  of  force  (or,  equivalently,  force  per  unit 
spanwise  length)  in  the  stream  wise  (or  _y-)  direction  for  the  control  volume  can  be 
written  as  [cf.  Eq.  (28)] 


where  rlv  is  the  wall  shear  stress  in  the  v-direetton  exerted  on  the  surface  of 
the  building  in  the  ->  (or,  vertical-)  direction.  The  first  term  in  Fq. (29)  is  the 
contribution  to  the  drag  force  resulting  from  the  pressure  discontinuity  across  the 
solid  obstacle  in  the  control  volume  and  the  second  term  is  the  contribution  arising 
from  the  frictional  force  along  the  surface  of  the  building!  s)  in  the  control  volume. 
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The  drag  force  FDx  enables  a  bulk  drag  coefficient  to  be  defined  as  follows: 


CD  = 


Ifbxj 

p  f;:0H  <u>2(=)  d= 


(30) 


where 


i  px=XQ+2H 

<«>(-)  =  2Hj  u(x,=)dx.  (3D 

Furthermore,  with  the  understanding  that  the  drag  force  always  acts  against  the 
direction  of  the  flow,  we  have  used  the  absolute  value  of  FDx  (i.e.,  magnitude  of  the 
drag  force  FDv)  in  the  definition  of  CD  to  ensure  that  CD  is  always  positive  definite. 
For  simplicity,  the  diagnosis  of  CD  (from  a  high-resolution  CFD)  using  the 
definition  in  Eq.  (30)  with  the  drag  unit  decomposition  for  the  obstacle  array  shown 
in  Fig.  10(a)  will  be  referred  to  henceforth  as  Method  (a). 

The  values  of  CD  for  each  drag  unit  in  Fig.  10(a)  obtained  using  Method  (a)  are 
shown  in  Fig.  11(a).  Here,  the  drag  coefficient  is  presented  in  the  form  of  a  non- 
dimensional  drag  parameter  CqAH /  V,  where  the  frontal  area  A  of  the  drag  unit  per 
unit  spanwise  length  is  A  =  FI  and  the  volume  V  of  the  drag  unit  per  unit  spanwise 
length  is  V  =  2 FI2,  so  A  =  AfV  =  (2 H)~l.  These  results  demonstrate  explicitly  that 
the  drag  coefficients  CD  of  the  first  two  drag  units  on  the  upstream  side  of  the  array 
are  significantly  larger  than  CD  of  the  remaining  four  drag  units  on  the  downstream 
side  of  the  array.  The  reduction  of  the  drag  coefficient  for  the  drag  units  “deep” 
within  the  array  may  be  an  example  of  the  ‘shelter  effect'  [41]  arising  from  the  mutual 
interference  of  the  canopy  elements. 

We  are  interested  also  in  the  determination  of  the  in-situ  drag  coefficient  Cd  of  a 
drag  unit  within  the  2D  obstacle  array  where  the  flow  becomes  fully  developed  (viz., 
in  the  region  of  the  array  where  there  is  negligible  streamwise  development  of  the 
flow  statistics).  In  order  to  determine  CD  for  the  streamwise  equilibrium  (fully 
developed)  condition  in  the  array,  we  extended  the  streamwise  extent  of  the  array  to 
include  14  2D  obstacles  (or,  13  drag  units).  Results  for  CD  from  the  high-resolution 
numerical  simulation  for  this  case  are  displayed  in  Fig.  11(b).  These  results  imply 
that  the  in  situ  drag  coefficient  CDAH  (or,  equivalently,  CDAH/V)  for  a  canopy 
element  within  the  array  in  the  regime  of  streamwise  equilibrium  is  about  1.67. 

The  major  disadvantage  of  the  drag  force  representation  exhibited  in  Fig.  10(a)  is 
that  a  drag  unit  is  associated  with  each  obstacle  in  the  array  and  the  drag  coefficient 
for  each  of  these  drag  units  is  required  as  input  information,  with  the  drag 
coefficients  determined  from  a  numerically  expensive  high-resolution  CFD  simula¬ 
tion.  The  computational  cost  increases  as  the  number  of  obstacles  in  the  array 
increases.  One  may  argue  that,  if  a  high-resolution  CFD  simulation  is  already 
available,  there  is  no  need  for  the  drag  force  approach.  Fortunately,  in  practice,  it  is 
not  necessary  to  associate  a  drag  unit  with  every  single  obstacle  in  the  array  because 
CD  approaches  its  fully  developed  value  within  the  array  relatively  quickly  as  shown 
in  Fig.  1 1(b),  although  CD  here  appears  to  deviate  slightly  from  its  fully  developed 
value  again  for  drag  units  12  and  13  on  the  downstream  end  of  the  array.  This  is  due 
to  the  pressure  effect  arising  from  the  wake  downstream  of  the  last  obstacle  (i.e.,  the 
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Fig.  II.  The  variation  of  the  non-dimensional  drag  parameter  CDAH/V  for  (a)  decomposition  of  the 
original  array  into  six  drag  units  and  (b)  decomposition  of  an  extended  array  into  thirteen  drag  units. 

14th  obstacle  of  the  array)  propagating  upstream.  The  results  of  Fig.  1 1(b)  provide 
the  motivation  for  the  development  of  an  alternative  method  for  the  representation 
of  2D  array  of  obstacles  by  the  drag  force  approach  as  illustrated  in  Fig.  10(b).  Here, 
the  first  and  last  obstacles  of  the  array  are  resolved  explicitly,  whereas  all  the 
intervening  obstacles  between  the  first  and  last  obstacles  are  represented  using  one 
drag  unit  oflength  \  \H  and  height  H.  The  drag  parameter  CDAH  for  this  single  unit 
is  taken  to  be  the  in  situ  drag  coefficient  in  the  regime  of  streamwise  equilibrium 
within  the  (extended)  array.  The  drag  representation  of  the  building  array  shown  in 
Fig.  10(b)  with  the  single  drag  unit  assigned  the  drag  parameter  CDAH  =  1  67  will 
be  referred  to  as  Method  (b)  henceforth. 
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4.3.  Results  and  discussion 

Because  the  representation  of  obstacles  using  a  distributed  drag  force  approach  in 
the  flow  model  “smears”  out  the  local  heterogeneity  of  the  detailed  flow  field  in  an 
obstacle  array,  it  is  convenient  (and,  perhaps,  more  appropriate)  to  compare  the 
predictions  of  mean  flow  and  turbulence  from  these  models  with  spatially  averaged 
flow  quantities.  As  an  example  of  these  flow  quantities,  Figs.  12(a)  and  (b)  show, 
respectively,  vertical  profiles  of  observed  mean  streamwise  velocity  u  and  turbulence 
kinetic  energy  k  at  10  x-positions  within  drag  unit  4  [cf.  Fig.  10(a)]  along  with  the 
horizontally  averaged  mean  velocity  and  turbulence  kinetic  energy  profiles  within  the 
drag  unit  (which  are  obtained  by  averaging  at  each  height  z,  u(x,z)  and  k{x,z)  over 
the  10  x-locations  within  the  drag  unit).  The  10  x-locations  shown  in  Fig.  12 
correspond  to  5  positions  over  the  rooftop  of  the  fourth  building  [xe(900, 1050)  mm] 
and  5  positions  over  the  fourth  street  canyon  [xe(1050, 1200)  mm]  of  the  array. 

The  local  heterogeneity  of  the  mean  flow  and  turbulence  energy  over  drag  unit  4  is 
clearly  seen  in  Fig.  12.  It  is  noted  that  the  horizontal  heterogeneity  of  k  is  stronger 
than  that  of  u.  In  particular,  profiles  of  k  below  the  building  height  at 
z  =  H  =  150  mm  near  the  lee  side  of  the  street  canyon  give  substantially  higher 
values  of  the  TKE  than  locations  within  the  protected  zone  near  the  windward  side 
of  the  canyon.  Unless  otherwise  indicated,  all  the  simulations  using  the  distributed 
drag  force  approach  shown  here  are  based  on  modifications  of  the  K-L  k-s 


Unit  4  Unit  4 


Fig.  12.  Profiles  of  the  mean  streamwise  velocity  fi  and  turbulence  kinetic  energy  k  measured  at  10  v- 
locations  within  drag  unit  4  (i.e.,  over  the  fourth  building  rooftop  and  the  fourth  street  canyon  of  the 
array).  The  horizontally  averaged  mean  streamwise  velocity  and  turbulence  kinetic  energy  over  the  drag 
unit  are  also  displayed. 
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turbulence  model  (by  the  inclusion  of  additional  source/sink  terms  in  the  k -  and  e- 
equations  as  described  in  Section  4.1). 

Fig.  13  exhibits  the  modeled  and  measured  mean  streamwise  velocity  and 
turbulence  kinetic  energy  obtained  with  the  distributed  drag  force  approach  at 
x  =  225  mm  (i.e.,  at  the  center  location  of  the  first  street  canyon).  At  this  location, 
note  that  the  horizontally  averaged  mean  streamwise  velocity  observed  experimen¬ 
tally  and  predicted  using  a  high-resolution  CFD  simulation  (denoted,  respectively,  as 
Expt-ave.  and  HR-ave.  in  Fig.  13)  are  in  very  good  agreement.  The  modeled  mean 
velocity  profiles  obtained  by  Method  (a)  differs  from  the  observed  values 
(represented  by  the  filled  circles  in  Fig.  13)  by  20%  or  less  within  the  roughness 
sublayer  (I  <z/H<2),  with  the  greatest  discrepancy  below  the  building  height 
z  =  H  =  150  mm  where  the  modeled  mean  flow  velocities  are  seen  to  under-predict 
the  magnitude  of  the  observed  mean  flow  reversal  near  the  ground.  Furthermore,  the 
modeled  mean  velocity  profiles  exhibit  the  greatest  shear  near  the  top  of  the  obstacle 
array,  with  greatly  reduced  shear  within  the  “urban”  canopy  which  accords  well  with 
the  observations.  Finally,  the  modeled  mean  velocity  profiles  obtained  using  Method 
(b)  (i.e.,  one  drag  unit)  provides  slightly  poorer  predictions  for  u. 

Because  the  turbulence  energy  levels  are  underestimated  by  the  high-resolution 
CFD  simulations,  it  is  not  surprising  that  the  horizontally  averaged  TKE  obtained 
using  the  high-resolution  CFD  also  underestimates  the  experimentally  observed 
horizontally  averaged  TKE  (cf.  Fig.  13).  Interestingly,  the  prediction  of  the 


U(m/s)  k  (m2/s2) 

Fig.  13.  Mean  streamwise  velocity  a  and  turbulence  kinetic  energy  k  profiles  predicted  using  the 
distributed  drag  force  approach  at  position  *  =  225  mm  (center  of  the  first  street  canyon)  are  compared 
with  observed  values  (Expt-ave.)  and  high -resolution  numerical  simulation  predictions  (HR-ave)  of  the 
horizontally  averaged  u  and  k. 
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horizontally  averaged  TKE  using  the  distributed  drag  force  approach  is  better  than 
that  obtained  by  using  the  high-resolution  CFD  simulation  (the  latter  of  which  is 
obtained  by  averaging  the  local  TKE  vertical  profiles  over  the  first  drag  unit).  Again, 
it  is  noted  that  the  modeled  TKE  profile  obtained  using  Method  (b)  provides  a 
slightly  poorer  prediction  than  that  obtained  using  Method  (a).  Nevertheless,  the 
simulated  TKE  for  both  methods  is  largest  within  the  roughness  sublayer  at 
z/Hx4/3,  and  is  reduced  significantly  within  the  “urban”  canopy  layer,  in  good 
agreement  with  the  observations. 

Figs.  14  and  15  compare  the  mean  streamwise  velocity  and  turbulence  kinetic 
energy  simulated  using  the  distributed  drag  force  approach  with  the  observed  values 
at  single  locations  x  =  525  mm  (center  of  the  second  street  canyon)  and 
x  =  1 125  mm  (center  of  the  fourth  street  canyon),  respectively.  The  mean  velocity 
profile  at  x  =  525  mm  for  Method  (a)  is  in  very  good  agreement  with  the 
experimental  results,  with  the  result  from  Method  (b)  giving  only  slightly  poorer 
agreement  in  this  case.  The  TKE  profiles  predicted  with  Methods  (a)  and  (b)  at 
x  =  525  mm  agree  quite  well  with  the  measured  profiles,  with  an  under-prediction  of 
the  turbulence  energy  levels  at  z/Hk 4/3  by  about  20%  at  most.  In  contrast,  at 
x  =  1 125  mm,  the  strong  shear  in  the  mean  velocity  profile  near  the  canopy  top  at 
z  —  H  is  under-predicted,  whereas  the  turbulence  energy  levels  within  the  urban 
canopy  (O^z/H ^\)  and  the  roughness  sublayer  ( 1  ^z/ II ^2)  are  over-predicted 
(the  more  so  within  the  urban  canopy  layer)  by  both  methods. 


Fig.  14.  Mean  streamwise  velocity  u  and  turbulence  kinetic  energy  k  profiles  predicted  using  the 
distributed  drag  force  approach  at  position  x  =  525  mm  (center  of  the  second  street  canyon)  are  compared 
with  observed  values  (Expt-ave.)  and  high-resolution  numerical  simulation  predictions  (HR-ave.)  of  the 
horizontally  averaged  u  and  k. 
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Fig.  15.  Mean  streamwise  velocity  u  and  turbulence  kinetic  energy  k  profiles  predicted  using  the 
distributed  drag  force  approach  at  position  *  =  1125  mm  (center  of  the  fourth  street  canyon)  are 
compared  with  observed  values  (Expt-ave.)  and  high-resolution  numerical  simulation  predictions  (HR- 
ave.)  of  the  horizontally  averaged  it  and  k. 


The  simulated  and  observed  mean  streamwise  velocity  and  turbulence  kinetic 
energy  at  x  =  1425  mm  (center  of  the  fifth  street  canyon)  are  compared  in  Fig.  16.  As 
for  the  case  at  x  =  1125  mm  (cf.  Fig.  1 5),  the  magnitude  of  the  mean  wind  shear  at  or 
near  z  =  H  is  under-predicted  by  both  methods.  However,  note  that  the  magnitude 
of  the  mean  velocity  at  positions  within  the  canopy  (z/H<  2/3)  is  in  good 
conformance  with  the  experimental  results.  Note  that  as  one  progresses  “deeper” 
into  the  obstacle  array,  the  differences  in  the  mean  velocity  predicted  by  Methods  (a) 
and  (b)  become  smaller  (as  the  mean  flow  converges  towards  the  condition  of 
streamwise  equilibrium  within  the  array).  Even  so,  the  condition  of  exact  streamwise 
equilibrium  for  the  spatially  averaged  time-mean  flow  within  the  canopy  has  not 
been  achieved  even  by  row  6  of  the  array.  In  particular,  note  from  Fig.  16  that  the 
spatially  averaged  time-mean  flow  is  slightly  negative  below  about  z/H<>  1  /4,  which 
agrees  with  the  reverse  flow  observed  here  in  the  spatially  averaged  experimental 
data  (Expt-ave)  and  predicted  by  the  spatially  averaged  high-resolution  CFD  results 
(HR-ave)  for  the  mean  velocity.  In  the  dense  2D  building  array  studied  here,  the 
vertical  momentum  flux  resulting  from  the  turbulent  stresses  do  not  penetrate  very 
deeply  into  the  urban  canopy  [in  particular,  note  from  Fig.  19  that  turbulent  shear 
stress  divergence  is  approximately  zero  near  the  base  of  2D  building  array 
(0<z///<;  1/2)],  so  the  flow  here  is  effectively  vertically  decoupled  from  the  flow 
aloft.  Deep  within  the  urban  canopy,  the  spatially  averaged  time-mean  momentum 
budget  reduces  to  a  balance  between  the  streamwise  advection  on  the  one  hand  and 
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Fig.  16.  Mean  streamwise  velocity  u  and  turbulence  kinetic  energy  k  profiles  predicted  using  the 
distributed  drag  force  approach  at  position  .*  =  1425  mm  (center  of  the  fifth  street  canyon)  are  compared 
with  observed  values  (Expt-ave.)  and  high-resolution  numerical  simulation  predictions  (HR-ave.)  of  the 
horizontally  averaged  «  and  k. 


the  pressure  gradient  and  canopy  drag  on  the  other  hand.  At  row  6  of  the  array,  the 
pressure  gradient  term  in  the  streamwise  spatially  averaged  momentum  transport 
equation  is  small,  but  not  exactly  zero  (implying  that  the  condition  of  streamwise 
equilibrium  here  is  not  exact).  This  adverse  pressure  gradient  force,  though  small 
within  the  canopy,  is  nevertheless  large  enough  to  cause  a  small  reverse  velocity  deep 
within  the  canopy  (z///<  1/4)  (viz.,  the  magnitude  of  the  spatially  averaged  time- 
mean  flow  tends  to  be  small  there,  so  even  a  very  weak  adverse  pressure  gradient 
force  can  result  in  a  back  flow  there). 

Note  that  the  spatially  averaged  time-mean  velocity  u  obtained  with  the 
distributed  drag  force  model  underpredicts  it  in  relation  to  both  the  HR-ave  and 
Expt-ave  (cf.  Fig.  16)  above  the  canopy.  In  addition,  the  shape  of  u  obtained  from 
the  distributed  drag  force  approach  is  different  than  that  exhibited  by  both  HR-ave 
and  Expt-ave  for  4/3<z///<7/3  (corresponding  roughly  to  the  inertial  sublayer). 
The  mean  flow  shown  in  Fig.  16  can  be  used  to  diagnose  a  roughness  length  z0  for 
2D  building  array,  which  in  this  case,  corresponds  to  a  skimming  flow  over  the  array. 
Fitting  u  in  the  inertial  sublayer  to  a  standard  semi-logarithmic  mean  velocity  profile, 
the  following  values  of  z0  were  inferred:  zo///*0.O085  for  Expt-ave;  zq/H *0.0095 
for  HR-ave;  and  z0/Hx 0.15  for  Method  (a)  and  Method  (b)  (which  are  obtained 
from  the  distributed  drag  force  approach).  Note  that  for  the  case  of  skimming  flow 
over  a  2D  array  of  buildings,  the  distributed  drag  force  approach  gave  a 
poor  estimate  for  zq.  However,  the  mean  flow  obtained  by  spatially  averaging  the 
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high-resolution  CFD  results  yields  a  value  for  z0  that  agrees  well  with  the  measured 
value  of  zo  (the  latter  obtained  from  the  spatially  averaged  experimental  data  for  the 
mean  velocity).  The  value  of  z<>  diagnosed  here  for  the  2D  building  array  from  the 
experimental  data  appears  to  be  larger  than  a  roughness  length  of  zo/Hx  0.001 
obtained  for  a  square  bar  array  (cf.  Fig.  2  in  [42])  with  a  roughness  density  of  A  —  0.5 
(defined  as  the  roughness  frontal  area  per  unit  ground  area).  This  difference  in 
roughness  length  may  be  the  result  of  the  roughening  of  2D  building  surfaces  (with 
grit-impregnated  paint)  in  the  wind  tunnel  experiment  of  [8],  in  order  to  prevent  re- 
laminarization  of  the  flow  on  the  building  surfaces.  For  skimming  flow  over  the  2D 
building  array  (analogous  to  a  square  bar  array),  where  the  normalized  displacement 
height  d/H >0.9  (which  reflects  the  fact  that  the  recirculation  cavity  between  the 
obstacles  is  isolated  from  the  main  flow  above  the  obstacle  array),  the  roughness 
length  z0  in  this  case  is  more  a  measure  of  the  momentum  absorption  by  the 

roughness  of  the  tops  of  the  obstacles,  rather  than  the  roughness  devolving  from 
the  individual  array  obstacle  elements  (square  bars)  themselves. 

At  x  =  1425  mm,  the  modeled  turbulence  energy  levels  from  both  methods  over¬ 
predict  the  measured  TKE  (by  a  factor  of  2  at  worst)  as  seen  in  Fig.  16.  The 
overestimation  of  TKE  in  the  region  of  the  obstacle  array  where  the  mean  flow 
within  the  canopy  has  become  fully  developed  (approximately  or  better)  may  be  due 
to  the  fact  that  the  additional  sink  term  in  Eq.  (25)  arising  from  the  presence  of  the 
obstacles  only  models  the  reduction  in  the  TKE  dissipation  arising  from  the 
reduction  in  the  TKE  as  the  turbulent  eddies  do  work  against  the  form  and  viscous 
drag  of  the  canopy  elements.  The  increase  in  the  dissipation  rate  due  to  the  reduction 
of  the  turbulent  mixing  length  scale  within  the  canopy  has  not  been  properly 
modeled.  It  is  conceivable  that  more  comprehensive  models  for  the  inclusion  of 
source/sink  terms  in  the  k-  and  e-equations  (as  advocated  by  Antohe  and  Lage  [24] 
and  Getachew  et  al.  [23])  may  lead  to  better  predictions  in  this  case  (albeit  at  the 
expense  of  much  more  complicated  models). 

Fig.  17  presents  the  distribution  of  the  turbulence  energy  levels  in  the  x—z  plane 
predicted  by  the  high-resolution  CFD  approach  and  the  distributed  drag  force 
approach  based  on  Methods  (a)  and  (b).  All  the  model  predictions  exhibit  a  zone 
with  large  production  of  TKE  centered  but  displaced  slightly  upwards  of  the 
windward  edge  of  the  roof  of  the  first  building.  However,  for  the  distributed  drag 
force  approach,  the  TKE  generated  in  this  zone  is  “exported”  (by  mean  advection 
or  turbulent  diffusion)  much  more  readily  into  the  “urban"  canopy  than  for 
the  high-resolution  CFD  simulation  where  the  obstacles  are  explicitly  resolved. 
Hence,  a  TKE  “plume”  extends  downstream  of  the  zone  of  maximum  TKE 
production,  but  the  damping  of  the  TKE  with  downstream  distance  and  within 
the  canopy  appears  to  be  too  small  for  the  two  distributed  drag  force  models. 
This  problem  is  most  serious  for  Method  (b)  where  all  the  obstacles  between  the 
first  and  last  buildings  are  represented  as  a  single  drag  unit.  However,  in  Method 
(a),  the  large  values  of  Co  assigned  to  the  first  and  second  drag  units  (cf.  Fig.  1 1) 
result  in  significant  dissipation  of  TKE  within  the  “urban”  canopy  with  the  result 
that  the  TKE  at  these  locations  are  actually  under-predicted  by  the  model  (cf.  Figs 
13  and  14). 
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Pig  1 7.  Turbulence  kinetic  energy  Lsopleths  over  the  2D  array,  obtained  using  a  high-resolution  numerical 
simulation  and  distributed  drag  force  approach,  are  compared 


4,4.  Dispersive  stresses  in  the  2D  building  array 

The  high-resolution  CFD  simulations  of  the  time-averaged  flow  within  and  over 
the  2D  building  array  (cf-  Section  3)  can  be  used  to  diagnose  the  dispersive  stresses  in 
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this  array.  The  dispersive  stresses  in  the  urban  canopy  arise  from  the  spatial 
correlations  in  the  point-to-point  variations  in  the  time-averaged  (mean)  velocity 
field.  The  dispersive  stresses  can  be  diagnosed  from  the  time-mean  velocity  field 
determined  from  the  high-resolution  CFD  simulation  as  follows.  In  this  section  only, 
it  is  necessary  to  distinguish  explicitly  the  time  averaging  operation  from  the  spatial- 
averaging  operation.  To  this  end,  overbars  will  be  used  to  denote  time  averaging 
whereas  angular  brackets  will  be  used  to  denote  spatial  averaging.  The  dispersive 
normal  and  shear  stresses  for  the  2D  building  array  [cf.  Fig.  10(a)]  can  be  diagnosed 
from  the  time-mean  velocity  field  w„  respectively,  as  follows: 

<  u"u"  >  (2)  =  1 5><a)(z)  -  <  u  >  (z))2  (32) 

a=l 

and 


1  N 

<  u"w"  > (z)  =  —  5>(«,(z)  -  <  u  > (z))(iv(a)(z)  -  <h->  (z)).  (33) 

a-[ 

In  Eqs.  (32)  and  (33),  u  =  u\  and  w  =  iij  are  the  streamwise  and  vertical  mean 
velocities,  respectively;  z  is  the  height  above  the  ground  surface  at  which  the  velocity 
is  evaluated;  a  double  prime  is  used  to  denote  a  departure  from  the  spatial  average; 
N  is  the  number  of  vertical  profiles  of  time-mean  velocity  in  a  drag  unit  [cf.  Fig. 
10(a)]  used  in  the  construction  of  the  spatial  average;  and  the  Greek  subscript  a 
enclosed  in  round  brackets  is  used  to  denote  a  running  index  for  the  vertical  profiles 
of  time-mean  velocity  in  a  drag  unit  (and,  should  not  be  confused  with  the  Roman 
subscript  index  /  used  to  distinguish  a  component  of  the  velocity).  The  angular 
brackets  in  Eqs.  (32)  and  (33)  refer  explicitly  to  a  spatial  average  over  a  drag  unit  of 
the  urban  array.  More  specifically,  to  resolve  the  variation  of  the  dispersive  and 
spatially  averaged  Reynolds  stresses  in  the  vertical  (or  z-)  direction,  the  averaging 
volume  implied  by  the  angular  brackets  is  taken  to  be  an  infinitesimally  thin 
horizontal  slab  through  the  drag  unit. 

Figs.  18  and  19  exhibit  the  normal  and  shear  dispersive  stresses,  respectively,  for 
drag  units  #1  to  #6  in  the  2D  building  array.  The  normal  and  shear  dispersive 
stresses  have  been  compared  with  the  spatially  averaged  normal  and  shear  Reynolds 
stresses  (denoted  by  <«V>  and  <n' wi,>,  respectively,  where  the  single  prime 
indicates  a  departure  from  the  time-averaged  value).  The  quantities  <^V>  and 
<nV>  represented  by  the  solid  lines  in  Figs.  18  and  19  were  determined  by  spatially 
averaging  the  Reynolds  stresses  (obtained  from  the  high-resolution  CFD  simulation) 
over  thin  horizontal  slabs  through  the  various  drag  units  (or,  averaging  units)  shown 
in  Fig.  10(a).  Interestingly,  for  the  2D  building  array,  the  normal  and  shear 
dispersive  stresses  are  comparable  in  magnitude  to  the  associated  spatially  averaged 
Reynolds  stresses  in  drag  unit  #  1  for  1  <  z/H  <  2,  but  appear  to  be  negligible  above 
canopy  height  ( z/H  >  1)  in  the  drag  units  further  downstream.  Furthermore,  the 
dispersive  shear  stress  here  is  opposite  in  sign  to  the  Reynolds  shear  stress.  For  drag 
units  #3  to  #  6,  the  normal  dispersive  stress  is  larger  than  the  normal  Reynolds  stress 
within  the  urban  canopy  for  the  height  ranges  2/3  £z/H<>  1  and0<z///<  1  /3.  Note 
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Fig.  18.  Comparison  between  the  Reynolds  norma!  stress  <wV>  (solid  line)  and  the  dispersive  normal 
stress  <u"u">  (dashed  line)  in  the  2D  building  array  for  drag  units  #1  to  #6. 


154 


F.S.  Lien  el  al  i  J.  Wind  Eng.  Ind.  Aerodyn.  92  (2004)  117-158 


CJomparison  be,ween  the  Reynolds  shear  stress  <  uV  >  (solid  line)  and  the  dispersive  shear  stress 
<“  &  >  (dashed  line)  in  the  building  array  for  drag  units  #1  to  #6. 
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that  for  the  2D  building  array  with  skimming  flow  over  the  array,  both  the  spatially 
averaged  Reynolds  normal  and  shear  stresses  tend  to  be  very  small  within  the  urban 
canopy  (z/H<\).  For  the  2D  building  array  in  the  skimming  flow  regime,  both  the 
dispersive  and  spatially  averaged  Reynolds  shear  stresses  are  negligible  within  the 
urban  canopy. 


5.  Summary  and  conclusions 

A  number  of  variants  of  the  k-e  turbulence-closure  model,  applied  with  wall 
functions,  has  been  used  to  calculate  the  disturbed  mean  flow  and  turbulence 
through  and  above  an  array  of  two-dimensional  buildings.  It  should  be  emphasized 
that  for  these  numerical  simulations,  the  closure  constants  in  the  k~e  turbulence 
model  retained  their  original  (standard)  values.  In  consequence,  all  model 
predictions  were  obtained  without  adjustment  of  any  of  the  closure  constants.  This 
ensures  that  the  process  did  not  simply  degenerate  into  a  curve-fitting  exercise. 
The  primary  purpose  of  this  effort  was  to  assess  critically  the  model’s  predictive 
performance  by  comparison  of  the  mean  flow  and  turbulence  statistics  with  a  very 
comprehensive  and  detailed  wind  tunnel  simulation.  In  addition,  the  results  of  the 
high  resolution  CFD  simulation  based  on  the  k~e  turbulence  models  were  used  also 
to  diagnose  drag  coefficients  that  can  be  used  in  a  distributed  drag  force 
representation  of  the  obstacle  array.  The  utility  of  using  a  distributed  drag  force 
representation  of  an  obstacle  or  group  of  obstacles  in  the  array  for  the  prediction  of 
the  disturbed  flow  has  been  investigated. 

The  comparison  of  predicted  (using  high-resolution  CFD  simulations)  and 
measured  streamwise  velocity  profiles  indicated  that  most  of  the  qualitative  features 
of  the  mean  flow  field— such  as  speed-up  over  the  first  building  rooftop,  the  presence 
of  a  thin  separation  layer  near  the  surface  of  the  first  rooftop  and  the  absence  of  such 
a  feature  in  all  the  succeeding  rooftops,  the  standing  vortex  in  the  street  canyon,  the 
thin  strong  shear  layer  that  forms  at  building  height,  the  recirculation  zone  behind 
the  last  building,  the  development  of  the  mixing  layer  above  and  downstream  of  this 
recirculation  zone,  and  the  rate  of  momentum  recovery  in  the  far  wake  region— are 
correctly  predicted.  The  quantitative  agreement  is  also  fairly  good  throughout  most 
of  the  flow  domain  for  the  mean  streamwise  velocity. 

The  turbulence  kinetic  energy  from  the  high-resolution  numerical  and  physical 
simulations  show  striking  resemblances  and,  indeed,  most  of  the  qualitative  features 
of  the  turbulence  energy  as  it  develops  through  and  above  the  array  of  buildings  are 
adequately  predicted.  In  terms  of  quantitative  agreement,  it  is  found  that  the  TKE 
level  is  generally  underestimated  by  the  model.  This  underestimation  may  be  caused 
by  the  inability  of  the  model  to  capture  the  disproportionate  and  inhomogeneous 
sensitivity  of  the  Reynolds  stress  components  to  secondary  strains  in  the  flow, 
particularly  those  associated  with  curvature.  Specifically,  the  production  of  TKE  just 
upstream  and  above  the  first  building  arising  from  the  mean  streamline  concave 
(destabilizing)  curvature  strain  here  and  from  the  “flapping”  of  the  shear  layer,  is 
not  adequately  captured  by  the  model.  Because  the  TKE  is  exported  downstream 
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and  vertically  by  advection  and  by  turbulent  and/or  pressure  diffusion,  the  under¬ 
prediction  of  the  TKE  at  and  near  the  leading  edge  of  the  first  building  (i.e.,  the 
location  where  the  turbulence  energy  production  is  greatest)  will  result  in  reduced 
TKE  levels  at  all  positions  downstream  of  this  location. 

It  has  been  shown  how  a  high-resolution  CFD  simulation  model  can  be  applied  to 
determine  the  drag  force  in  a  given  control  volume  of  the  obstacle  array  and  how  this 
information  can  be  used  to  diagnose  the  values  of  the  drag  coefficient  Cd  required  in 
a  distributed  drag  force  representation  of  the  obstacles  in  the  array.  These  diagnosed 
values  of  CD  have  been  used  in  various  drag  force  models  of  the  obstacle  array,  and 
have  been  compared  to  detailed  measurements  of  the  horizontally  averaged  mean 
velocity  and  turbulence  energy  levels  in  the  array.  It  was  found  that  the  numerical 
simulations  based  on  the  distributed  drag  force  approach  were  able  to  satisfactorily 
predict  most  of  the  main  characteristics  of  the  mean  streamwise  velocity,  and  many 
of  the  relevant  features  of  the  TKE.  However,  although  many  of  the  qualitative 
features  of  the  flow  around  and  within  the  obstacle  array  were  correctly  simulated 
with  the  distributed  drag  force  approach,  the  complicated  details  relating  to  the  local 
heterogeneity  of  the  flow  cannot  be  predicted.  For  the  latter,  a  high-resolution  CFD 
simulation  is  required. 

The  standard  or  K-L  k-e  turbulence-closure  model  with  an  isotropic,  linear  eddy 
viscosity  is  perhaps  the  simplest  complete  turbulence  model  (viz.,  no  advance 
knowledge  of  any  property  of  the  turbulence  is  required  for  the  simulation  other 
than  the  initial  and/or  boundary  conditions  for  the  problem)  that  is  available 
currently,  and  its  moderately  good  predictive  performance  of  the  complicated 
developing  flow  through  an  array  of  two-dimensional  buildings  without  the  need  to 
adjust  any  closure  constants  is  encouraging.  This  model  may  be  useful  as  a  general 
purpose  simulator  of  urban  flows  since  it  is  simple  enough  to  be  tractable 
numerically  and,  hence,  not  require  excessive  computing  time.  However,  if  a  better 
representation  of  mean  velocity  and  turbulence  fields  is  required  in  the  vicinity  of  a 
cluster  of  buildings,  then  the  non-linear  k—e  model  is  recommended.  The  latter 
model  is  slightly  more  complicated  than  the  standard  k — e  model,  and  requires 
between  about  10  and  20  percent  more  computational  effort  than  the  standard  k—e 
model  to  provide  a  prediction  for  the  flow.  It  is  conceivable  that  these  numerical 
simulation  models  could  provide  all  the  statistics  of  the  disturbed  wind  flow  in  a 
building  cluster  required  as  input  to  a  physically  based  model  for  the  prediction  of 
dispersion  of  contaminants  in  the  array. 
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